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An  analysis  of  cicada  mating  calls,  measured  in  field  experiments,  indicates  that  the  very  high  levels 
of  acoustic  energy  radiated  by  this  relatively  small  insect  are  mainly  attributed  to  the  nonlinear 
characteristics  of  the  signal.  The  cicada  emits  one  of  the  loudest  sounds  in  all  of  the  insect 
population  with  a  sound  production  system  occupying  a  physical  space  typically  less  than  3  cc.  The 
sounds  made  by  tymbals  are  amplified  by  the  hollow  abdomen,  functioning  as  a  tuned  resonator,  but 
models  of  the  signal  based  solely  on  linear  techniques  do  not  fully  account  for  a  sound  radiation 
capability  that  is  so  disproportionate  to  the  insect’s  size.  The  nonlinear  behavior  of  the  cicada  signal 
is  demonstrated  by  combining  the  mutual  information  and  surrogate  data  techniques;  the  results 
obtained  indicate  decorrelation  when  the  phase-randomized  and  non-phase-randomized  data 
separate.  The  Volterra  expansion  technique  is  used  to  fit  the  nonlinearity  in  the  insect’s  call.  The 
second-order  Volterra  estimate  provides  further  evidence  that  the  cicada  mating  calls  are  dominated 
by  nonlinear  characteristics  and  also  suggests  that  the  medium  contributes  to  the  cicada’s  efficient 
sound  propagation.  Application  of  the  same  principles  has  the  potential  to  improve  radiated  sound 
levels  for  sonar  applications.  [DOI:  10.1 121/1.3050258] 

PACS  number(s):  43.60.Wy,  4 3.25. Ts,  43.80. Jz  [EJS]  Pages:  958-967 


I.  INTRODUCTION 

A.  Background 

The  objective  of  this  research  is  to  begin  to  understand 
how  the  cicada,  a  small  insect,  emits  one  of  the  loudest 
sounds  in  all  of  the  insect  population  despite  its  relatively 
small  size.  Detailed  knowledge  of  the  characteristics  of  the 
cicada's  acoustic  signature  is  a  necessary  step  toward  the 
ultimate  goal  of  transferring  this  biotechnological  feat  of  na¬ 
ture  to  a  manmade  transduction  system  of  similar  propor¬ 
tions.  The  cicada’s  highly  effective  sound  production  system 
occupies  a  physical  space  typically  less  than  3  cc.  Cicadas 
are  sexually  dimorphic,  and  only  males  possess  the  structures 
necessary  for  making  loud  audible  sounds.  Male  sounds  are 
broadcast  advertisements  for  attracting  females,  and  they  are 
typically  loud,  rhythmic,  and  easily  distinguished  from  back¬ 
ground  noise.  Males  create  sound  by  flexing  a  pair  of  ridged 
abdominal  membranes  called  tymbals.  The  sounds  made  by 
these  tymbals  are  amplified  by  the  hollow  abdomen  function¬ 
ing  as  a  tuned  resonator,  as  described  by  current  research. 
Nevertheless,  the  tuned  resonator  explanation  in  the  current 
literature  does  not  account  for  the  sound  radiation  capabili¬ 
ties  of  the  cicada. 

Studying  the  sound  production  system  of  the  cicada  in 
captivity  has  some  inherent  difficulties:  cicadas  vocalize  only 
on  sunny  days  and  do  not  respond  as  well  to  indoor  lighting. 
However,  the  recent  discovery  of  a  female  response  to  the 
mating  call  of  the  male  has  made  it  possible  to  conduct  ex¬ 
periments  on  cicadas  in  the  field.  In  many  species,  the  fe- 
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males  answer  loud  male  signals  with  quiet  wing  flick 
responses.  Males  perceiving  such  responses  will  ap¬ 
proach  that  signal  and  continue  to  call  even  if  disturbed 
These  female  signals  are  easily  imitated,  which  provides  an 
important  tool  for  collecting  and  manipulating  cicadas:  in  an 
acoustical  duet  with  a  female,  a  male  will  become  sexually 
excited  and  continue  to  sing  even  if  disturbed  or  manipu 
lated.  Similar  manipulations  without  a  duet  may  cause  a 
male  to  stop  calling. 

B.  Current  experimental  opportunities  with  periodic 
cicadas 

Opportunities  to  collect  and  study  cicadas  of  the  mid- 
western  and  eastern  United  States  are  limited  by  the  insects’ 
periodical  cycles.  The  insects  have  a  13  or  17  year  life  cycle 
and  emerge  in  mass  numbers,  known  as  broods,  at  predict¬ 
able  times  in  predictable  locations.  Since  the  sounds  and 
behaviors  of  the  cicadas  are  not  as  yet  fully  characterized, 
an  upcoming  emergence  of  the  insects  will  provide  an  oppor¬ 
tunity  to  test  and  report  on  the  tymbal  and  abdomen  struc¬ 
tural  dynamics  that  generate  the  cicada’s  mating  call.  Mea¬ 
suring  the  cicada’s  two  most  important  anatomic  structures 
for  producing  sound  would  add  to  the  scientific  understand¬ 
ing  of  the  extent  of  the  nonlinear  nature  of  the  cicada  signals 
and  would  also  help  to  explain  the  high  sound  levels  pro¬ 
duced  by  this  small  insect. 

C.  Current  state  of  cicada  research 

A  comprehensive  review  of  the  prominent  journal  ar¬ 
ticles  on  cicada  sound  production  and  mechanisms  and  a 
survey  of  a  number  of  subject  area  expert  textbooks  in  this 
field  revealed  that  the  mechanisms  underlying  the  cica¬ 
da's  sound  levels  and  efficient  sound  propagation  are  not 
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fully  understood.  The  cicada  song  has  been  classically  mod¬ 
eled  using  linear  mathematical  methods.  However,  these  lin¬ 
ear  methods  are  insufficient  for  a  true  model  of  the  system 
because  the  buckling  tymbals  within  the  cicada  sound  pro¬ 
duction  system  are  essential  to  the  acoustic  level  and  propa¬ 
gation  of  this  mating  call.  Inelastic  buckling  is  well  recog¬ 
nized  as  a  nonlinear  phenomenon  among  researchers. 

II.  QUANTIFYING  NONLINEARITY  IN  CICADA  SIGNALS 

A.  Technical  approach 


nal  signal  x(k).  Nonlinearity  exists  if  the  randomized  signal’s 
Mi  diverges  from  the  original  signal’s  Mi,  thereby  signifying 
that  nonlinearity  must  be  present  in  the  time  series. 

Linear  calculations  are  plotted  with  the  nonlinear  signal 
results  in  order  to  display  the  contrast  between  linearity  and 
nonlinearity.  The  Gaussian  rule  in  Eq.  (3)  is  presented  to 
indicate  how  the  correlation  coefficient  and  Mi  arc  linked  for 
the  special  case  of  a  pair  of  Gaussian  random  variables  with 
normalized  correlation  coefficient  p: 

Ia.B-~2  1°S2(1  -P2)-  (3) 


A  signal  processing  repertoire  includes  methods  to  test 
for  (a)  Gaussianity,  (b)  non-Gaussianity,  (c)  linearity,  and  (d) 
nonlinearity.  The  basic  analytical  tools  used  to  perform  these 
tests  arc  the  temporal  power  spectrum,  the  average  mutual 
information  (Mi)  (i.e.,  an  information  theory  technique),  and 
surrogate  data  hypothesis  testing  (i.e.,  phase  randomization). 
A  preliminary  determination  of  the  degree  and  influence  of 
these  effects  in  the  sound  production  system  of  the  cicada  is 
explored  using  the  Volterra  expansion. 


B.  Nonlinear  signal  processing  with  mutual 
information 


In  order  to  substantiate  the  existence  of  nonlinearity,  a 
quantitative  method  must  be  established.  This  study  used  Mi 
(Ref.  30)  and  a  surrogate  data  method  to  confirm  nonlinear¬ 
ity.  Equation  ( 1 ),  which  defines  the  general  Mi,  is  a  probabi¬ 
listic  equation  used  quantitatively  to  assess  information  be¬ 
tween  two  random  variables  A  and  B. 


'a.B=X 

Q 


log2 


(i) 


where  PAB(ai^Y  is  the  joint  probability  of  events  from  sets 
A={a,}  and  /?={/?,},  and  PA{a^  and  PB(bj)  are  the  marginal 
individual  probabilities  associated  with  sets  A  and  B,  respec¬ 
tively.  For  example,  set  B  can  be  taken  as  a  collection  of 
events  that  are  time-delayed  versions  of  the  events  in  set  A, 
which  is  the  case  for  this  research. 

The  surrogate  data  method  consists  of  randomizing  the 
phase  of  a  signal  spectrum  as  shown  in 

N-\ 

s(n)  =  2  x(k)c~l2m,klN  for  n  =  0 :N  -  1 . 

1=0 

(2) 

Af/2-1 

S(k)  =  X  s(n)ei['p(")+2’rat/A'] 

n=0 

N-l 

+  x  s(n)e^n)+2mk,N^  for  k=0:N-  1 . 

n=N!2 


S(k)  is  the  overall  phase-randomized  discrete  Fourier  trans¬ 
form  (DFT)  of  signal  x(k),  where  phases  {^>(/i)}  are  indepen¬ 
dent,  uniform,  and  randomly  distributed  over  a  2 tt  range. 
s(/i)  is  the  complex  amplitude  spectrum  of  the  DFT  of  the 
original  time  series,  which  is  altered  by  a  random  phase 
tp{n).  In  the  surrogate  method,  the  inverse  DFT,  S(&),  is  cal¬ 
culated  from  s(/i),  and  this  transformed  time  series  is  called 
the  surrogate  data  and  used  for  a  comparison  with  the  origi¬ 


Note  that  the  Mi  1A  B  is  0  when  the  correlation  coefficient  p  is 
0,  while  the  Mi  goes  to  infinity  as  the  correlation  coefficient 
goes  to  positive  or  negative  1 .  This  holds  for  all  joint  Gauss¬ 
ian  processes,  which  are  considered  linear  processes.  The 
frequency  domain  equivalent  of  the  correlation  coefficient 
for  measuring  linearity  is  the  subject  of  many  published  pa¬ 
pers  on  coherence  in  the  1993  reprint  text  by  Carter.  Papers 
on  coherence  include  how  to  estimate  coherence,  how  to 
minimize  bias  and  variance,  and  how  to  determine  confi¬ 
dence  bounds  for  estimates  of  magnitude-squared  coherence. 
In  general,  for  stationary  random  process,  proper  averaging 
of  large  time  segment  improves  estimation. 


C.  Higher-order  spectral  techniques  using  magnitude- 
squared  bicoherence 


Furthermore,  multispectral  techniques  exist,  such  as  bi¬ 
coherence  and  tricoherence,  which  could  provide  additional 
understanding  of  a  non-Gaussian  process.  If  the  cicada  sig¬ 
nals  are  non-Gaussian,  the  bicohcrence  could  determine  if  a 
process  is  a  mixed-phased  process  or  a  nonlinear  process. 
Cumulants  are  higher-order  statistical  information  about  any 
data  series.  For  example,  the  first-order  cumulant  for  a  sta¬ 
tionary  process  is  its  mean  value.  The  second-order  cumulant 
for  a  zero-lag  process  is  the  covariance,  the  third  is  skew¬ 
ness,  and  the  fourth  is  kurtosis.  As  the  first  step  beyond  first- 
order  spectral  analysis,  this  study  analyzes  the  bicoherence. 
Equation  (4)  is  the  magnitude-squared  bicoherence  (MSB), 
which  is  used  to  provide  evidence  on  whether  a  data  se¬ 
quence  is  linear  or  nonlinear: 


\S2(fuf2)\2 


(4) 


D.  Nonlinear  signal  processing  using  a  Volterra 
expansion 

A  nonlinear  fit  is  performed  on  the  cicada  data  by  means 
of  the  Volterra  expansion,  which  describes  the  first-order  and 
second-order  signal  dynamics  present  within  the  cicada  time 
series.  The  data  acquired  on  the  cicada  consist  of  two  laser 
measurements,  denoted  by  sequences  {a (n)}  and  {b(n)|,  as 
well  as  a  simultaneously  sampled  microphone  sequence 
{z(/i)J.  This  situation  will  be  considered  to  be  a  two-input, 
onc-output,  nonlinear  “cicada  system”  with  mathematical 
memory,  instead  of  the  traditional  system  of  one  input  and 
one  output.  This  formulation  lends  itself  to  a  Volterra  ex- 
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pansion,  which  can  be  used  to  determine,  quantitatively,  the 
extent  of  nonlinearity  present  between  the  two  inputs  and  the 
one  output. 

1.  Second-order  Volterra  formulation 

The  standard  Volterra  expansion  for  a  onc-input  system 
takes  the  form 

*i-i 

y(rt)  =  h0+  X  h{(k)a(n-k) 

4=0 

*2-l  *2-1 

+  S  S  h2(kj)a(n  -  k)a(n  -/),  (5) 

4=0  7  =4 

when  carried  to  the  second  order,  where  y(n)  is  the  Volterrra 
fit.  The  three  functions  h0,h {(k),  and  h2(kj)  are  the  zeroth- 
ordcr,  first-order,  and  second-order  kernels,  respectively.  The 
first-order  terms  are  carried  out  to  length  Kly  while  the 
sccond-ordcr  terms  are  carried  out  to  length  K2  by  K2.  Be¬ 
cause  of  the  symmetry  inherent  to  h2  in  Eq.  (5),  the  summa¬ 
tion  index  j  can  be  limited  to  value  k  and  above.  The  un 
known  kernels  appear  linearly  in  model  Eq.  (5),  whereas  the 
known  excitation  {a (n)}  appears  nonlinearly  through  a  prod¬ 
uct  of  delayed  versions. 

With  a  two-input  model,  a  generalization  is  necessary, 
namely, 

y(rt)  =  yo+yi('0  +  y2('IK  (6) 

where  components 

*1-1  *1-1 

yi(n)=  X  hB(k)a(n  —  k)+  X  hb(*)b(/i  -  A), 

4=0  4=0 

(7) 

* 2-l  *2-l 

y2('0=  X  X  h^(kj)a(n -k)a(n-j) 

4=0  7=4 

*2"1  *2"! 

+  I  I  hbb(kj)h(n  -  k)b(n  -  j) 

4=0  y=4 

*2-l  *2-l 

+  X  2  hab(*,y)a(/j-A:)b(«-y). 

4=0  7=0 

There  are  two  linear  components  in  {y,(/j)}  each  of  length 
AT, ,  and  three  nonlinear  (second-order)  components  in  model 
output  {y2(/;)}.  The  advantage  of  the  inherent  symmetry  in 
the  two  auto  components  in  y2(n)  reduces  the  number  of 
kernel  values  that  have  to  be  determined.  However,  the  cross 
component  h ab(Jt,y)  in  y2(n)  has  no  such  symmetry  and 
therefore  requires  a  full  K2  by  K2  expansion.  The  total  num¬ 
ber  of  unknown  kernel  coefficients  in  Eqs.  (6)  and  (7)  is 

K=  1  +  2AT,  +  K2(2K2  +  1).  (8) 

It  is  desired  to  choose  these  coefficients  h  so  that  the  total 
model  output,  Eq.  (6),  fits  the  measured  microphone  output 
z{n)  as  well  as  possible  using  least  squares,  so  that  y(/i) 
~z (n).  The  least-squares  approach  is  adopted  because  the 
simultaneous  equations  for  the  optimum  kernel  coefficients 
will  then  all  be  linear. 


Although  the  laser  and  microphone  data  have  been 
sampled  at  frequency  fs= 96  kHz,  the  crucial  frequency  con¬ 
tent  of  the  cicada  system  itself  is  not  believed  to  extend 
above  10  kHz.  Therefore,  reductions  in  the  memory  lengths 
and  K2  for  the  first-order  and  second-order  kernels  in  the 
model  will  not  alias  the  cicada  statistical  information  Also, 
the  memory  lengths  are  decreased  to  minimize  the  computer 
random  access  memory  required  to  calculate  Eq.  (61.  Conse¬ 
quently,  a  decimation  factor  of  M  is  applied  to  the  kernels. 
Thus,  the  model  to  be  fitted  is  a  modification  of  Eq.  (7), 
namely, 

*H  *i-> 

y  i(/i)  =  X  h*(k)a(n  -  M k)  +  X  hb(k)b(n-  Mk) 

4=0  4=0 


=  ya('l)  +  yb(«), 

*2-1  *2-1 

y2(n)  =  X  X  hu(£,./)a(fl  -  Mk)a(n  -  Mj) 

4=0  7=4 
* 2-l  A 2-l 

+  E  S  hbb(*,y)b(n  -  Mk)b(n  -  Mj) 

4=0  7=4 
*2-l  *2-l 

+  2  X  babU,y)a(/i  -  Mk)b(n  -  Mj) 

4=0  7=0 

“  y«»(«)  +  ybb(«)  +  y.b  ("). 


(9) 


By  this  means,  the  memory  length  of  the  first-order  kernels  is 
MK\  units  of  the  sampling  increment  1  lfs,  while  that  of  the 
second-order  kernels  is  MK2  units.  Notice  that  the  measured 
data  a (n),  b(/i ),  and  z (n)  are  not  decimated,  thereby  retaining 
any  harmonic  and  intermodulation  products  that  might  have 
been  created  by  the  cicada  system  itself. 


2.  Least-squares  considerations 

The  details  of  a  first-order  fitting  procedure  will  be  pre¬ 
sented;  this  formulation  can  then  be  extended  to  include  all 
the  terms  in  Eq.  (9).  The  pertinent  equation  that  governs  the 
least-squares  approach  is  to  make 

*1-1 

ya(")=  2  h,(A)a(n  -  Mk) 

4=0 

=  ha(0)a(/i)  +  ha(  1  )a (n  -  M) 

+  •••  +ha(tf,  -  \  )a(n  -  M(K}-  1))  (10) 

approximate  z(/i)  for 

/V,  <  /i  <  Nx  =  M(K\  -  1)  +  1 ,  (11) 

where  N,  is  the  common  data  length  of  the  three  available 
data  sequences.  The  particular  starting  value  A,  for  n  arises 
so  that  the  inherent  buildup  transient  of  the  first-order  kernel 
ha(/c)  will  be  excluded  from  the  fitting  procedure  using  Eq. 
(10).  Trying  to  fit  the  transient  can  only  degrade  the  proce¬ 
dure;  confining  the  error  minimization  to  the  steady-state 
model  output  is  the  best  approach. 

Equations  (7)  and  (9)  can  be  put  into  a  matrix  formula¬ 
tion  as  follows: 
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a(N )) 

a(Ni-M)  • 

••  a(l) 

a(Nl+l) 

a(  2) 

a(N,) 

••  a(N,-Nl  +  \)_ 

MO) 

z(N ,) 

MO 

~ 

z(Ni  +  1) 

_M*i-0_ 

_z(Nt) 

In  matrix  notation,  this  equation  reads 

Dh  ~  Z,  (13) 

where  matrix  D  is  (Nt-Nx  +  1 )  X  column  vector  h  is 
X  1,  and  column  vector  Z  is  (N,-N{  +  1 )  X  1  Since  data 
length  Nt  is  generally  a  very  large  number,  whereas  the  num¬ 
ber  A'j  of  kernel  coefficients  is  usually  small,  the  attempted 
fit  in  Eqs.  (12)  and  (13)  cannot  be  achieved  exactly.  If  an 
error  sequence  is  defined  as  the  difference  between  the  left- 
hand  and  right-hand  sides  of  Eq.  (13),  and  the  sum  of 
squared  errors  is  minimized,  it  can  be  shown  that  the  opti¬ 
mum  kernel  h0  satisfies  the  equations 

D'Dh0  =  D'Z,  h0  =  (D'D)-'D'Z  =  (D'D)  \  (D'Z). 

(14) 

Observe  that  both  sides  of  Eq.  (13)  are  premultiplied  by  D' 
and  the  approximation  is  replaced  by  an  equals  sign. 

In  MATLAB  notation,  the  least-squares  solution  to  Eq. 
(13)  is  obtained  according  to 

h0  =  D  \  Z.  (15) 

The  approach  in  Eq.  ( 15)  instead  of  Eq.  (14)  is  more  advan¬ 
tageous,  in  that  the  condition  number  of  square  matrix  D'D 
is  the  square  of  the  condition  number  of  /),  which  makes 
approach  Eq.  ( 14)  less  reliable. 

111.  CICADA  EXPERIMENTS 

A.  Cicada  field  experiment  with  microphone 

Measured  data  from  the  Tibicen  chloromera  and  Tibicen 
lyriccn  species  were  recorded  in  the  field  with  a  parabolic 
microphone  and  a  Marantz  PMD670  behind  a  high  school  in 
West  Hartford,  CT.  Thus,  some  noise  from  crickets  and  other 
environmental  sounds  contaminated  the  signals,  but  this 
noise  was  about  40-60  dB  below  the  cicada’s  distinct  call. 
Also,  only  distinctive,  strong  cicada  calls — not  low-level 
idling  songs — were  used  to  determine  the  existence  of  non¬ 
linear  behavior  in  the  signal. 

B.  Cicada  field  experiment  with  laser  and  microphone 

Another  field  experiment  was  conducted  on  Magicicada 
septendecim  and  Magicicada  cassini  in  Julibce  State  College 
Park  in  Peoria,  IL  by  Hughes  and  Katz  as  well  as  the  ac¬ 
knowledged  support.  The  test  site  was  specifically  chosen  to 
amass  live  specimens  during  the  17  year  emergent  cycle  of 
periodical  cicadas  in  brood  XIII.  In  this  test,  the  tymbal  mo¬ 
tion  of  the  cicada  was  measured  by  laser  while  the  acoustic 
output  was  recorded  via  a  microphone.  This  simultaneous 
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measurement  was  made  with  dual  Polytec  OFV-508  optical 
measurement  heads  controlled  by  the  Poly  tec  electronic  sig¬ 
nal  processor  OFV-2802.  The  dual  optical  sensor  heads  al¬ 
lowed  simultaneous  measurement  of  the  motion  of  both  tym¬ 
bals  and  the  tymbal-to-abdomen  motion.  Meanwhile,  the 
parabolic  microphone  was  used  to  continuously  record  either 
the  output  of  the  two  tymbals  or  a  tymbal-abdomen  experi¬ 
mental  setup.  The  measurements  provide  the  quantitative  ve¬ 
locity  from  the  laser,  which  can  be  transformed  into  displace¬ 
ment  (position  versus  time)  information.  The  tymbal  motion 
is  proportional  to  acoustic  mechanical  vibration  and  thus  can 
be  used  in  conjunction  with  the  microphone  acoustic  output 
data.  Thus,  the  data  gathered  in  this  field  experiment,  in 
which  the  output  and  input  signals  were  obtained  simulta¬ 
neously. ,  allows  a  unique  opportunity  to  gain  further  scientific 
understanding  of  the  cicada  sound  production  system.  The 
input  signal  is  from  the  tymbal,  and  the  output  signal  is  re¬ 
corded  via  the  microphone.  These  data  sets  allow  the  appli¬ 
cation  of  nonlinear  mathematical  techniques,  such  as  the  Vol- 
terra  expansion,  to  real-world  signals. 

Two  channels  were  designated  for  collecting  dual  laser 
measurements  of  the  motion  of  both  tymbals  and  the  tymbal- 
abdomen  simultaneously.  Most  collections  consisted  of  live 
wingless  insects,  otherwise  intact.  A  third  channel  was  used 
to  record  the  microphone  output.  A  fourth  channel  was  used 
to  time  tag  all  signals.  Two  multichannel  recorders  were  used 
to  record  and  back  up  the  data. 

The  experimental  setup  used  a  multichannel  digital  re¬ 
corder  set  to  a  96  kHz  sampling  rate,  a  laser  and  its  control¬ 
ler,  a  backup  digital  recorder,  and  a  digital  timekeeper.  The 
digital  recorder  facilitated  the  simultaneous  acquisition  of  the 
tymbal  motion,  a  microphone,  and  time  stamp.  One  advan¬ 
tage  afforded  by  this  experimental  arrangement  was  the  op¬ 
portunity  to  analyze  quantitatively  both  synchronous  and 
asynchronous  tymbal  motion  and  how  such  motion  might 
impact  the  vocalization  output 

IV.  NONLINEAR  SIGNAL  PROCESSING  RESULTS 

A.  Nonlinear  and  linear  modeling  of  cicada  signal 
results 

A  linear  representation  of  a  simulated  time  series  of  the 
cicada  vocalization  based  on  the  superposition  of  a  Gaussian 
white-noise  signal  passed  through  three  parallel  independent 
narrowband  filters  and  the  power  spectral  representation  of 
this  simulated  time  series  representation  is  shown  in  Fig. 
1(a).  Figure  Kb)  shows  the  temporal  power  spectrum  of  the 
cicada  vocalization  measured  in  the  field.  Given  the  striking 
similarity  between  these  two  spectra,  the  immediate,  but  er¬ 
roneous,  conclusion  might  be  that  the  modeled  spectra  are  a 
good  representation  of  the  insect’s  actual  acoustic  signature. 
This  is  the  classical  error  made  by  analysts  using  strictly 
linear  techniques  to  model  behavior  in  a  physical  system.  In 
fact,  linear  techniques  are  insufficient  for  a  correct  analysis 
of  the  underlying  signal  processing  mechanics  of  the  cicada’s 
sound  production  system. 

The  simulated  signal’s  Mi  (i.e.,  white  Gaussian  excita¬ 
tion  passed  through  three  parallel  narrowband  filters)  and  the 
Gaussian  Mi  model  of  the  simulated  signal  [from  Eq.  (  3)]  are 
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FIG.  1 .  (Color  online)  Power  spectra  of  simulation  vs  field  measurement,  (a) 
Power  spectrum  of  simulated  cicada  call  estimated  with  linear  model,  (b) 
Temporal  power  spectrum  of  cicada  call  as  measured  in  the  field  The  spec¬ 
tra  appear  similar,  hut  the  linear  model  does  not  accurately  represent  the 
measured  acoustic  signature 


FIG.  3.  (Color  online)  The  Gaussian  curve  fit  to  the  prohahility  density 
function  (dashed  line)  is  very  similar  to  the  simulated  data  curve  hased  on  its 
histogram  (solid  line) — textbook  example  of  results  for  a  linear  Gaussian 
process  and  baseline  for  the  field  measurements  of  cicada  signals. 


compared,  which  produce  an  identical  plot  to  Fig.  2  where 
the  solid  line  is  the  Mi  and  the  “X”  or  cross  symbol  indicates 
the  Gaussian  rule.  However,  Fig.  2  compares  a  plot  of  the  Mi 
obtained  from  the  non-phase-randomized  simulated  date 
with  the  phase-randomized  simulated  data  As  expected,  the 
phase-randomized  signal  and  the  non-phase-randomized  sig¬ 
nal  Mi  generate  the  same  values  since  phase  randomization 
implies  near  Gaussianity  (according  to  the  central  limit  theo¬ 
rem),  which  in  turn  implies  linearity.  In  Fig.  3,  the  Gaussian 
curve  fit  to  the  probability  density  function  (PDF)  is  very 
similar  to  the  simulated  data  curve  based  on  its  histogram. 
These  graphs  arc  like  textbook  examples  of  what  the  results 
should  be  for  a  linear  Gaussian  process.  The  graphs  also 
serve  as  a  baseline  with  which  to  compare  the  field  measure¬ 
ments  of  cicada  signals.  A  field  recording  of  the  cicada  signal 
is  displayed  in  Fig.  4(a),  with  a  zoomed-in  version  of  the 


FIG.  2.  (Color  online)  Mi  obtained  from  the  phase-randomized  simulated 
signal  as  a  function  of  the  time  lag  r  (X)  compared  with  the  non-phase- 
randomized  simulated  data  (solid  line  plot)  The  marker  (X)  is  indistinguish¬ 
able  from  the  solid  line,  as  expected,  since  phase  randomization  implies  near 
Gaussianity  (according  to  the  central  limit  theorem),  which  in  turn  implies 
linearity.  If  plotted,  the  simulated  signal’s  Mi — white  Gaussian  excitation 
passed  through  three  narrowband  filters— as  a  function  of  the  time  lag  r 
(solid  line)  is  compared  with  the  Gaussian  rule  model  of  the  simulated 
signal  from  Eq.  (3)  (x).  The  plots  are  identical  to  what  is  shown. 
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same  signal  in  Fig.  4(b).  The  expanded  view  illustrates  the 
complexity  of  the  cicada  call.  In  Fig.  5,  the  cicada  data  (as 
Mi  of  the  cicada  vocalization)  are  compared  with  the  surro¬ 
gate  data  (the  phase-randomized  cicada  data)  to  show  how 
the  phase-randomized  cicada  signal  approaches  the  Gaussian 
rule.  The  separation  between  Mi  plots  of  the  cicada  signal 
and  the  phase-randomized  version  of  the  same  signal  cor 
roborate  its  non-Gaussian  probability  distribution.  Also,  the 
separation  between  the  Mi  of  the  cicada  vocalization  and  the 
simulated  data  using  the  Gaussian  rule  is  identical  to  Fig.  5, 
which  is  consistent  with  strong  non-Gaussianity  [i.e.,  Eq. 

m 

Additional  evidence  of  the  non-Gaussian  behavior  of  the 
cicada  vocalization  is  illustrated  in  Fig.  6.  This  figure  com¬ 
pares  the  Gaussian  fit  with  the  curve  based  on  the  histogram 
associated  with  the  cicada's  non-Gaussian  signal.  Note  that 
the  non-Gaussian  behavior  in  the  cicada  call  was  not  re¬ 
flected  in  the  modeled  spectra  shown  earlier  in  Fig.  1  Thus, 
Fig.  6  reinforces  the  point  made  earlier  against  relying  solely 
on  power  spectral  techniques  to  analyze  cicada  time  wave¬ 
forms. 
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FtG.  4.  (Color  online)  Time  scries  of  measured  cicada  vocalization:  a  0.5  s 
field  recording  (a)  and  an  expanded  segment  illustrating  the  complexity  of 
the  cicada  call  (b). 
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FIG.  5.  (Color  online)  Mi  of  the  cicada  vocalization  as  a  function  of  the 
time  lag  r  (solid  line)  compared  with  the  surrogate  data,  i.e,.  Mi  of  the 
phase-randomized  cicada  data  (X).  The  separation  between  Mi  plots  of  the 
cicada  signal  and  the  phase-randomized  version  of  the  same  signal  confirm 
its  non-Gaussian  probability  distribution.  If  plotted,  the  cicada  signal’s  Mi  as 
a  function  of  the  time  lag  r  (solid  line)  is  compared  with  the  Gaussian  rule 
model  of  the  signal  from  Fq.  (3).  (X)  is  identical  to  this  plot. 


B.  Results  of  higher-order  spectral  technique  using 
MSB 

The  results  from  Eq.  (4)  are  shown  in  this  section  of  the 
report.  Namely,  the  bicoherence  indicates  whether  the  time 
series  data  is  linear  or  nonlinear.  The  theoretical  MSB 
B(f\,f2)  of  a  linear  process  is  flat  throughout  the  /j-/2 
plane.  Therefore,  the  example  shown  in  Fig.  7  with  an  exact 
solution  is  utilized  to  ascertain  if  the  bicoherence  will  always 
predict  a  flat  surface  for  this  linear  process,  where  the  center 
frequency  ( fc )  is  zero.  Here  is  a  zero-mean,  unit-variance 
Gaussian  process.  Consider  the  entire  expression  of  Eq.  (4) 
with  the  denominator  terms  containing  the  individual  fre¬ 
quencies  as  well  as  the  addition  of  both  frequencies.  Figure  8 
illustrates  that  the  denominator  terms  do  not  completely  flat¬ 
ten  the  conical  peak  when  the  numerator  (bispectrum)  is  di¬ 
vided  by  the  denominator  shown  in  Eq.  (4).  Note  that  the 
MSB  plot  is  not  flat  suggesting  that  the  process  is  not  linear. 
Figure  9  illustrates  that  the  PDF  does  have  an  effect  on  the 
random  amplitudes  generated  in  the  bicoherence  computa- 


FIG.  6.  (Color  online)  The  PDF  of  the  cicada  signal  (solid  line,  multiple 
peaks)  compared  with  a  Gaussian  curve  fit  to  the  PDF  (dashed  line,  single 
peak).  The  separation  of  the  plots  is  additional  evidence  of  the  non-Gaussian 
behavior  of  the  cicada  signal. 


FIG.  7.  (Color  online)  The  spectrum  of  nonlinear  function  x(t)2-  1  provides 
an  opportunity  to  assess  the  capability  of  the  MSB  to  evaluate  nonlinearity. 


tion.  However,  the  amplitudes  of  the  sample  MSB  in  Fig.  9 
also  suggest  that  averaging  may  help  reduce  the  peaks  and 
valleys  of  the  plots.  Figure  10  confirms  that  averaging  does 
reduce  the  amplitude  variation  of  the  sample  MSB.  Such 
observations  are  consistent  with  the  bias  and  variance  reduc¬ 
tion  observed  in  coherence  estimation.  Nonetheless,  the 
amplitude  variation  is  not  completely  removed  from  the  es¬ 
timate.  Hence,  averaging  the  MSB  only  scales  down  the 
problem  of  a  fluctuating  surface.  Consequently,  a  flatness 
factor  could  yield  a  meaningful  parameter  to  employ  on  the 
variations  calculated  for  the  sample  MSB. 
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Equation  (16)  defines  the  “flatness  factor,”  which  derives 
from  the  sample  standard  deviation  of  the  fluctuations  in  the 
f\  ~fi  plane.  The  variable  a  represents  the  calculation  of  the 
sample  standard  deviation  for  a  known  sample  set  {a,}  with  a 
total  sample  size  /V,  while  the  parameter  x  is  the  sample 
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FIG.  8.  The  amplitude  of  the  MSB  for  #(f)2-l  is  large  despite  the  denomi¬ 
nator  attempting  to  scale  the  numerator  (bispectrum) 
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FIG.  9.  The  MSB  is  calculaied  for  a  normalized  random  noise  signal  and  the 
mesh  ploi  exposes  that  surface  undulations  are  present  in  the  f\-fi  plane, 
which  would  give  the  false  conclusion  that  Gaussian  random  noise  is  not  a 
linear  process. 


FIG.  11  (Color  online)  A  CDF  for  an  exponential  random  noise  averaged 
over  1000  trials  indicates  the  Gaussian  characteristics  of  the  MSB. 


mean  of  the  data  set,  Therefore,  the  sample  standard  devia¬ 
tion  <tmsb  represents  the  variation  for  the  entire  plane. 

Sinee  the  PDF  of  a  given  data  set  affects  the  standard 
derivation,  a  determination  must  be  made  to  show  how  the 
PDF  affects  the  MSB.  A  cumulative  distribution  function 
(CDF)  is  a  convenient  approach  to  determine  if  a  PDF  has 
Gaussian  characteristics,  which  is  the  case  for  the  exponen¬ 
tial  and  normal  distributions.  If  the  exponential  and  the  nor¬ 
mal  CDFs  differ,  the  implication  is  that  the  PDF  influences 
the  MSB  significantly.  In  Fig.  1 1,  the  CDF  for  the  exponen¬ 
tial  signal  is  shown.  Compared  to  Fig.  12,  there  is  a  differ¬ 
ence  in  value  of  the  accumulation  of  bins  in  the  MSB  for  the 
exponential  and  normal  distributions  of  0.5-0.95,  respec¬ 
tively.  This  difference  is  considerable  and  thus  indicates  that 
the  PDF  of  the  data  contributes  notably  to  the  MSB.  Conse¬ 
quently,  a  prerequisite  algorithm  must  account  for  the  effect 
of  the  PDF  on  the  MSB,  and  PDFs  could  vary  between  in¬ 
dividual  cicadas  and  even  more  so  with  different  species. 
Because  the  MSB  significantly  depends  on  the  PDF,  addi¬ 
tional  work  beyond  the  scope  of  the  present  research  is  re¬ 
quired  for  quantifying  how  statistically  useful  the  sample 
MSB  is  for  analyzing  cicada  signals. ' 


C.  Volterra  expansion  method  applied  to  cicada 
experimental  results 

A  sample  of  100,000  data  points  from  the  laser  and  mi¬ 
crophone  measurements  is  selected  for  fitting  purposes,  us¬ 
ing  a  decimation  factor  Af=4  The  number  of  coefficients 
employed  for  the  first-order  fit  is  ^  =  100,  and  the  number 
used  for  the  second  order  is  AS  =  50,  which  results  in  solving 
for  a  total  of  K- 5251  coefficients.  The  execution  time  re¬ 
quired  simply  to  fill  the  D'D  matrix,  using  segment  length 
A=1000,  is  1670  s  on  a  2.4  GHz  computer.  The  solution 
time  for  the  optimum  kernel  is  22  s,  and  the  time  required  to 
compute  all  five  individual  component  waveforms 
y*(n),yb(n),yM(n),yhh(n),  and  yab(n)  in  Eq.  (4)  is  220  s. 
The  singular  value  decomposition  of  matrix  D'D  took  720  s, 
and  its  condition  number  is  7.2  X  106.  Thus,  approximately  8 
decimal  digits  (out  of  15)  of  significance  remain  in  the  nu¬ 
merical  results  obtained.  The  total  execution  time  is  almost 
44  min.  The  ratio  of  the  power  in  the  total  model  output,  per 
Eq.  (6),  to  the  power  in  the  measured  microphone  output 
z(n)  is  0.41.  Thus,  the  Volterra  fitting  procedure  of  the  sec¬ 
ond  order  is  capable  of  representing  41%  of  the  power  of  the 
microphone  waveform. 


FIG.  10.  The  mesh  plot  of  MSB  for  a  normalized  random  noise  averaged 
over  1000  trials  slill  has  ripples  on  lhe  /,-/2  plane.  Flowever,  this  averaging 
scales  the  fluctuating  surface  but  determining  nonlinearity  remains  difficult 
lo  quanlify. 


FIG.  12  (Color  online)  CDF  for  a  normal  random  noise  averaged  over  I  (XX) 
trials  is  also  Gaussian  but  has  a  significant  differenl  PDF,  which  affecls  lhe 
MSB  when  compared  to  the  exponential  random  noise  in  Fig.  1 1 . 
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FIG  13.  (Color  online)  Short  time  segments  of  simultaneous  laser  measure¬ 
ments  of  the  motion  of  the  two  tymbals  and  the  sound  output  at  the  micro¬ 
phone  in  the  field  are  shown.  The  repeated  bursts  of  narrowband  energy  are 
clear  in  the  laser  A  and  B  measurements  in  (a)  and  (b),  respectively.  The 
microphone  output  (c)  tends  to  capture  environmental  noise,  which  blurs  the 
individual  bursts  seen  in  (a)  and  (b). 


FIG  15.  (Color  online)  Least-squares  fit  for  Volterra  expansion,  (a)  A  seg¬ 
ment  of  the  microphone  data  (solid  line)  [see  Fig.  1 3(c)]  is  compared  with 
the  total  fitted  waveform  (dashed  line).  The  fit  is  rather  good  in  some  time 
intervals,  but  poorer  in  others — a  manifestation  of  the  fitted  power  ratio  of 
0.41.  (b)  and  (c)  are  the  individual  total  first-  and  second-order  fits,  respec¬ 
tively. 


As  a  check  case,  an  unrelated  random  white-noise  pro¬ 
cess  replaces  the  microphone  output  z{n ),  and  the  fitting  pro¬ 
cedure  is  repeated  with  identical  parameters.  The  power  ratio 
of  the  final  fit  is  reduced  to  0.053.  An  additional  run  with  a 
different  random  sequence  for  z (n)  yields  a  comparable 
value  for  the  power  ratio.  Thus,  the  fitting  ratio  0.41  that  is 
actually  attained  is  a  significant  value  and  indicates  that  non- 
linearities  are  present  in  the  cicada  system  between  laser  in¬ 
puts  and  microphone  output.  In  fact,  the  power  in  the 
second-order  component  y2(«)  in  Eq.  (0)  is  almost  three 
times  greater  than  the  power  in  the  linear  component  y\(n). 
Also,  the  power  in  the  cross  component  yab(/i)  is  greater 
than  the  powers  in  the  two  autocomponents  y aa(n)  and 
ybb("). 

D.  Volterra  graphical  results 

A  short  time  segment  of  the  three  simultaneous  field 
measurements  is  displayed  in  Fig  13.  Each  of  the  lasers — A 


Frequency  (Hz) 


FIG  14.  (Color  online)  Unsmoothed  spectral  estimates  of  the  data  measured 
at  lasers  A  (a)  and  B  (b)  and  at  the  microphone  (c).  These  complete  data 
segments,  corresponding  to  both  tymbal  and  microphone  outputs  in  Fig.  1  3, 
are  used  in  the  Volterra  fit. 


and  B — measured  the  narrowband  buckling  of  a  single  ci¬ 
cada  tymbal,  shown  in  Figs.  13(a)  and  13(b),  respectively, 
while  the  microphone  [see  Fig.  13(c)]  tended  to  capture  all 
the  acoustic  noise  in  that  segment  of  the  field  measurement. 
The  corresponding  (unsmoothed)  spectral  estimates  of  the 
complete  data  segments  used  in  the  Volterra  fit  are  plotted  in 
Fig.  14.  A  segment  of  the  microphone  data  zin)  is  compared 
with  the  total  fitted  waveform  y(n)  in  Fig.  15(a).  The  fit  is 
rather  good  in  some  time  intervals  but  poorer  in  others — a 
manifestation  of  the  actual  fitted  power  ratio  of  0.41.  The 
plots  in  Figs.  15(b)  and  15(c)  show  the  individual  total  first- 
order  and  total  second-order  fits,  respectively.  The  two  first- 
order  ttme-domain  kernels  are  plotted  in  Fig.  16,  while  their 
complex  frequency-domain  transfer  functions  are  displayed 
in  Fig.  17.  The  decimation  factor  M- 4  is  the  reason  for  the 
upper  frequency  limit  of  12  kHz  for  these  kernels. 

An  alterative  Volterra  expansion  is  computed  that  mea¬ 
sures  the  nonlinearity  in  the  cicada’s  mating  call.  Table  I 
describes  each  tymbal’s  contribution  to  the  first-and-second- 
order  components  in  the  cicada  signal  by  using  a  Volterra 
expansion  with  M  =  3  and  K2=l 0.  The  Volterra  technique 


FIG  16.  (Color  online)  First-order  time-domain  kernels  of  tymbal  motion 
measured  by  lasers  A  (a)  and  B  (b). 
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requires  both  output  and  input  functions:  in  this  case,  micro¬ 
phone  data  zz  are  the  output,  and  the  input  functions  are  laser 
A  and  laser  B.  Variable  y  in  Table  1  represents  the  Volterra 
expansion  solution.  The  lowercase  a  and  b  for  the  Volterra 
expansion  correspond  to  the  measured  vibration  from  the 
two  tymbals.  The  first-order  and  second-order  Volterra  solu¬ 
tions  are  y,  and  y2,  respectively.  Therefore,  ya  is  the  first- 
order  solution  for  laser  A,  yb  is  the  first-order  solution  for 
laser  B,  and  yaa  represents  the  second-order  solution  for  laser 
A,  and  so  forth.  The  combination  of  both  tymbals  yab  is  a 
mixed  second-order  component  that  contributes  approxi¬ 
mately  36%  of  the  second-order  Volterra  solution.  This  find¬ 
ing  suggests  that  parametric  generation  may  contribute  to  the 
cicada's  sound  propagation.  Also,  the  second-order  solution 
contains  87%  of  the  Volterra  solution,  which  also  indicates 
significant  nonlinearity  in  the  cicada  call. 

There  are  two  reasons  that  the  fitting  fraction  is  not 
greater  than  0.41  or  0.36  in  the  cases  above.  The  common 
length  Kx-  100  of  the  two  first-order  kernels  ha(A)  and  hb(A) 
is  apparently  adequate  because  both  of  the  estimated  first- 
order  kernels  decayed  essentially  to  zero  at  both  ends  of  the 
memory  interval  of  length  MKX.  Because  the  microphone 
data  contained  a  considerable  amount  of  background  noise, 
the  common  length  /G=50  of  the  three  second-order  kernels 
haa (kj),  hbb (kj),  and  hab(fc,y)  is  inadequate  because  the 
estimate  does  not  decay  sufficiently  by  the  edges  of  the  k 
-  j  plane.  The  common  length  =  did  reduce  the  decay  in 
the  kernel,  but  still  more  delay  is  required. 

E.  implications  of  cicada  signal  results  for  underwater 
acoustics 

In  the  commercial  world,  small  active  “fish-finding”  so¬ 
nars  with  wristband  receivers  are  sold  at  a  modest  price  to 
sports  fishermen.  Recent  scientific  advances  in  the  study  of 
small  insects  generating  loud  acoustics  in  air,  by 
Bcnnet-Clark,  suggest  opportunities  to  transfer  this  tech¬ 
nology  to  small-sized  active  sonar  applications.  The  cicada’s 
tymbal  mechanics  include  a  tymbal  plate,  four  ribs  and  rcsi- 
lin  pad.  “The  four  ribs  of  the  Timbal  buckle  inwards  from 
posterior  to  anterior....  A  train  of  four  sound  pulses  [was 
produced],  each  corresponding  to  the  inward  buckling  on  one 
rib....”  The  difference  in  the  frequencies  “suggests  that  the 
mass-to-stiffness  ratio  that  determines  the  resonant  frequen¬ 
cies  of  the  various  pulses  differs  from  pulse  to  pulse.... 
Pulses  produced  later  in  the  inward  buckling  sequence  were 
less  affected  by  the  loading  than  earlier  ones.  This  suggests 
that  the  effective  mass  determining  the  resonance  in  the  later 
pulses  is  greater  than  that  in  the  earlier  pulses. ...  The  Tymbal 
appears  to  act  as  an  energy  storage  mechanism  that  releases 
energy  as  the  tymbal  ribs  buckle  inwards  in  sequence....  It 
[also]  appears  that  the  central  part  of  each  rib  is  decoupled 
from  its  predecessor  as  it  buckles  inwards....  observations 
suggest  that  the  vibration  is  initially  non-linear  when  the 
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FIG.  17.  (Color  online)  Transfer  functions  of  firsi-order  lime-domain  ker¬ 
nels  of  lymbal  motion:  power  (a)  and  phase  (b)  for  laser  A  dala,  and  power 
(c)  and  phase  (d)  for  laser  B  dala 


amplitude  of  vibration  is  largest  but  that  it  becomes  more 
linear  as  the  pulse  decays.”  In  certain  regions  of  interest, 
Bennet-Clark  suggests  that  the  resonance  “is  determined  by 
the  simple  interaction  between  linear  mass,  compliance  and 
damping  elements”*  Moreover,  Young  and  Bennet-Clark 
suggest  “that  the  inward  movement  of  the  tymbal  plate  and 
rib  1  will  be  constrained  by  the  resilin  sheet  that  couple  it  to 
the  next,  unbuckled,  rib.”  Bennet-Clark  states  that  the 
“change  in  relative  resonance  frequency  implies  the  ratio  be¬ 
tween  the  mass  and  the  stiffness  of  the  vibrating  system  al 
ters  as  successive  ribs  buckle.”  Young  and  Bennet-Clark  go 
on  to  observe  that  “analogous  mechanisms  for  producing 
long  coherent  waveforms  by  a  succession  of  impulses  have 
been  modeled  in  bush  crickets.”  Part  of  the  cicada's  sound 
generation  relies  on  the  rubber-like  resilin  material  which 
“acts  as  an  energy  store  in  which  muscle  work  is  stored 
comparatively  slowly  (perhaps  during  the  first  2  to  4  ms  of 
contraction)  and  is  then  released  rapidly  by  the  sudden  buck¬ 
ling  of  the  tymbal  ribs  over  a  period  of  one  cycle  or 
230  microsec.”  Note  that  230  (is  is  0.23  ms,  and,  thus,  for 
2.3  ms  (between  2  and  4  ms),  the  discharge  occurs  in  nomi¬ 
nally  ten  times  less  time.  This  detailed  description  of  the 
cidada’s  storage  and  release  mechanism  provides  useful  in 
sight  for  construction  of  an  underwater  acoustic  source. 

V.  CONCLUSIONS 

This  research  utilizing  advanced  signal  processing  tech¬ 
niques  to  support  that  the  cicada’s  loud  mating  call  is  pro¬ 
duced  as  a  non-Gaussian,  nonlinear  vocalization.  Calculating 
the  Mi  to  determine  nonlinearity  is  often  difficult  because  of 
the  many  degrees  of  freedom  (i.e.,  large  sample  size)  re¬ 
quired  to  obtain  a  valid  joint  probability  function  [Eq.  (I )]. 
To  circumvent  size  limitations  on  the  data  sample,  the  Volt¬ 
erra  expansion  was  used  to  discover  that  the  second-order 
kernel  possesses  enough  of  the  energy  in  the  cicada  signal  to 


TABLE  I  Volierra  expansion  to  quantify  cicada  signal  nonlinearily  for  laser  A  (a),  laser  B  (b),  and  the  microphone  (zz). 


Ya 

Yb 

Yaa 

ybb 

yab  Y1 

Y2 

Y  zz  Error 

6.4677  X  10'K 

3.2090  x  10-7 

2.5461  X  1(T7 

1.0602  x  10-6 

6.0081  x  10"7  3.9690  x  10“7 

1 .6650  x  10* 

1 .9100  X  10  *  4.81 80  X  I0"6  2.9090  x  10“* 
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indicate  nonlinear  behavior.  Although  the  curse  of  dimen¬ 
sionality  pertains  to  the  Volterra  expansion,  the  current 
second-order  solution  is  satisfactory  for  initiating  the  quan¬ 
tification  of  the  nonlinear  parameters  for  modeling  and  simu¬ 
lating  purposes. 

Mi  and  surrogate  hypothesis  testing  techniques  can  be 
combined  with  the  Volterra  or  Wiener  expansion  to  param¬ 
etrize  the  cicada’s  abdominal  cavity  motion  and  tymbal  ex¬ 
citation  and  create  a  model  to  simulate  the  insect’s  sound 
production  mechanism.  The  combined  and  individual  contri¬ 
butions  of  the  motion  of  the  abdomen  and  tymbals  can  be 
evaluated  and  quantified  using  the  nonlinear  least-squares 
technique  combined  with  a  Wiener  expansion.  The  Wiener 
method  allows  the  computation  of  second-  and  third-order 
solutions  with  fewer  coefficients.  Reducing  the  coefficient 
size  for  each  order  makes  the  computational  limitation  less 
of  a  concern,  and  consideration  can  be  given  to  calculating 
higher-order  kernels.  Determining  these  parameters  will  aid 
in  developing  a  device  that  generates  sound  propagation  with 
the  efficiency  of  the  cicada  vocalizations.  Such  efficient 
sound  wave  propagation,  if  viable  in  water,  would  markedly 
enhance  source  radiation  efficiency  for  a  variety  of  sonar 
applications. 
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